1. Field of the Invention
The invention relates to radiation therapy planning for the treatment of tumors or for stereotactic radiosurgery and more particularly to the optimization of the radiation dose delivered to a patient by inverse treatment planning.
2. Description of the Background Art
The delivery of radiation for the treatment of tumors or obliteration of other abnormalities in human patients has undergone important changes in the last few years. The most important change has been the use of fine pencil beams of radiation of fixed intensity that can be scanned in a manner similar to a raster, remaining in each position of the scan for different lengths of time. With groups of pencil beams directed to the tumor or abnormality from different angles around the body, it is possible to deliver a radiation dose distribution that is effective in the treatment of the tumor or abnormality while sparing organs or tissues that are in the path of the radiation beams. That technology is known as Fluence Modulated Radiation Treatment (FMRT) or Fluence Modulated Radiation Surgery (FMRS). For historical reasons, it is also commonly known as Intensity Modulated Radiation Treatment/Surgery (IMRT/S) although the latter terminology is misleading since it implies modulated beam intensities.
The problem of therapy planning for FMRT/S is a difficult one: the fluences of a large number of pencil beams from many angles have to be calculated so that the sum of all their effects results in a uniform dose to the volume that contains the tumor or abnormality. Excessive dose to some portions of that volume can result in medical complications and under-dosing a region can lead to poor tumor control or recurrence. At the same time, the calculated fluences must fulfill some restrictive requirements to the radiation dose that is allowable in the different sensitive tissues traversed by the beams. The requirements placed by a therapist may be medically desirable but, to a smaller or larger extent, will be inconsistent with the physical laws that govern the absorption of radiation by the different tissues in the patient. The calculation of beam fluences, then, requires a methodology that allows for those inconsistencies and still delivers a set of beam flux results that achieve the most important of the desired results, as defined by the therapist. The process of obtaining that set of results is usually termed xe2x80x9coptimizationxe2x80x9d. The volume that includes the tumor or other abnormalities to be treated is normally termed xe2x80x9cPlanning Treatment Volumexe2x80x9d, or PTV. There can be more than one PTV in a treatment plan, but, without loss of generality, this document will refer to a PTV in the singular. Similarly, the volume that includes sensitive tissues or organs is normally termed xe2x80x9cOrgan at Riskxe2x80x9d, or OAR. There can be more than one OAR, but the term will be used in the singular, without loss of generality, except when giving specific examples with multiple OAR""s.
Optimization is invariably carried out by maximizing or minimizing a target or cost function incorporating the requirements and limitations placed by the therapist on the solution of the problem. Once a target or cost function has been defined, an algorithm is called upon to find the values of the beam fluence variables that will maximize or minimize that function, as the case may require.
Current art in optimization for FMRT/S can be divided into technologies that use:
a) a mathematically derived iterative formula to maximize or minimize the target or cost function. Those methods need an analytic target or cost function, so that it has explicit first partial derivatives with respect to the beam fluence variables.
b) a stochastic method that repeatedly tries maximizing or minimizing a target or cost function by testing whether a small change in the fluence of a randomly selected beam results in a change of that function in the desired direction.
Current art for the two types of technologies can be summarized as:
a1) using a Dynamically Penalized Likelihood target function which maximizes the likelihood that the resulting beam fluence values will yield the desired dose to the PTV, subject to a maximum dose specified for the OAR (see U.S. Pat. No. 5,602,892, J. Llacer, xe2x80x9cMethod for Optimization of Radiation Therapy Planningxe2x80x9d, Feb. 11, 1997; and J. Llacer, xe2x80x9cInverse radiation treatment planning using the Dynamically Penalized Likelihood methodxe2x80x9d, Med. Phys., 24, (11) pp 1751-1764, 1997).
a2) using a quadratic cost function which is to be minimized, with additional quadratic penalty functions imposed to apply restrictions. A gradient method speeded up by scaling the gradient with the inverse of the diagonal elements of the Hessian matrix is used for that minimization. It allows for specifying the maximum dose to the OAR and also the fraction of the OAR volume that is allowed to receive more than a certain dose (see T. Bortfeld, J. Bxc3xcrkelbach, R. Boesecke and W. Schlegel, xe2x80x9cMethods of image reconstruction from projections applied to conformation radiotherapyxe2x80x9d, Pys. Med. Biol., 1990, Vol 35, No. 10, 1423-1434; T. Bortfeld, J. Bxc3xcrkelbach and W. Schlegel, xe2x80x9cThree-dimensional solution to the inverse problem in conformation radiotherapyxe2x80x9d, Advanced Radiation Therapy Tumor Response Monitoring and Treatment Planning, Breit Ed., Springer, pp 503-508, 1992; T. Bortfeld, W. Schlegel and B. Rhein, xe2x80x9cDecomposition of pencil beam kernels for fast dose calculations in three-dimensional treatment planningxe2x80x9d, Med. Phys. 20 (2), Pt. 1, pp 311-318, 1993; and T. Bortfeld, J. Stein and K. Preiser, xe2x80x9cClinically relevant intensity modulation optimization using physical criteriaxe2x80x9d, Proceedings of the XIIth International Conference on the Use of Computers in Radiation Therapy (ICCR), Leavitt and Starkschall, eds., Springer, pp 1-4, 1997).
a3) using a quadratic cost function with constraints that limit the space of the allowable solutions to those that are non-negative (there cannot be negative beam fluences) and lead to doses to the OAR that are not above a maximum and/or result in no more than a specified fraction of the OAR volume receiving more than a certain dose. A modified form of the Conjugate Gradient method is used for minimization of the cost function (see S. V. Spirou and C-S. Chui, xe2x80x9cA gradient inverse planning algorithm with dose-volume constraintsxe2x80x9d, Med. Phys. 25 (3), pp 321-333, 1998).
b1) using the Simulated Annealing technique, a stochastic method, to minimize a variety of proposed cost functions, including some possibly important biological functions like Tumor Control Probability and Normal Tissue Complication Probability (see S. Webb, xe2x80x9cOptimizing radiation therapy inverse treatment planning using the simulated annealing techniquexe2x80x9d, Int. Journal of Imaging Systems and Tech., Vol. 6, pp 71-79, 1995, which summarizes the extensive work over many years by that author and co-workers) and
b2) using non-analytic cost functions that describe the fractions of PTV and OAR volumes that are to receive no more or less than a certain range of doses. The minimization of those cost functions is carried out by the simulated annealing method, (see U.S. Pat. No. 6,038,283, M. P. Carol, R. C. Campbell, B. Curran, R. W. Huber and R. V. Nash, xe2x80x9cPlanning method and apparatus for radiation dosimetryxe2x80x9d, Mar. 14, 2000).
The above methods are either in use at several institutions in the U.S. and elsewhere or being prepared for inclusion in Radiation Therapy Planning commercial software packages by several Companies.
An advantage attributed to the stochastic method of U.S. Pat. No. 6,038,283 is that it allows the therapist to define cost functions that prescribe the fractions of PTV and OAR volumes that are to receive no more or less than a certain range of doses. In order to provide that advantage, that method has to use non-analytic cost functions without explicit first partial derivatives. Testing whether a small change in the fluence of a randomly selected beam will result in an increase or decrease of the cost function involves complex calculations that slow down the optimization considerably. The number of times that those complex calculations have to be carried out for a prostate cancer case, for example, may be 5xc3x97105.
Short optimization calculation times are important in the clinic. Even with substantial experience, a therapist may have to recalculate a therapy plan several times, each time with a somewhat different set of PTV and OAR restrictions, until a satisfactory compromise optimization has been found.
Because of the frequent inconsistencies imposed by Physics between the doses that a therapist would want to prescribe for a PTV and an OAR, it is important not to over-constrain the problem. What is needed is an implementation of the Simulated Annealing method that allows the therapist to specify the fractions of the OAR volume that are to receive no more than a certain range of doses but allows the optimization algorithm to do its best for the PTV automatically. Furthermore, if the method uses an analytic cost function with partial derivatives that can be calculated rapidly, the therapist will have the ability to explore several possible optimizations within a clinically practical computation time and select the one that best suits the medical requirements of the patient.
The present invention is a new Adaptive Simulated Annealing (ASA) method of optimization. It uses the stochastic method of Simulated Annealing, in conjunction with dynamic characteristics related to the Dynamically Penalized Likelihood method of U.S. Pat. No. 5,602,892, to minimize an adaptive analytic cost function, with desired values of dose for each of the voxels of the OAR selected by a novel methodology. It allows the therapist to specify the fraction of OAR volumes that are to receive no more than a certain range of doses, while yielding a dose in each of the PTV voxels that is as close as physically possible to the desired dose.
Specifically the invention includes a method and apparatus for solving a numerical optimization problem that yields pencil beam fluences that will result in the optimum treatment plan using a predetermined set of pencil beams.
The preferred method and apparatus for solving the numerical optimization problem comprises a computer running a new ASA iterative algorithm. Taking into consideration the desired dose in all the PTV voxels, the range of maximum doses desired for different fractions of the OAR volume and the energy deposited per unit fluence into the patient""s tissues by all the selected pencil beams, the new ASA algorithm finds the beam fluences that minimize the cost function:       B    ⁢          (      a      )        =                    ∑                  i          ∈          D                    ⁢                        w          i                ⁢                  Φ          ⁢                      (                                          h                i                                  (                  k                  )                                            ,                              d                i                                      )                                +                  ∑                              i            ∈            S                                              (                                                h                  i                                      (                    k                    )                                                  -                                  s                  i                                            )                         greater than             0                              ⁢                        β          i                ⁢                  Θ          ⁢                      (                                          h                i                                  (                  k                  )                                            ,                              s                i                                      )                              
wherein
a=vector of pencil beam fluences
hi(k) is the dose received by voxel i at iteration k,
di=desired dose in voxel i of the PTV,
D=the region that includes all the voxels of the PTV,
si=desired dose in voxel i of the OAR,
S=the region that includes all the voxels of the OAR,
"PHgr"( )=an analytic function of the dose received by voxel i in the PTV at iteration k and of the desired dose for the same voxel,
"THgr"( )=an analytic function of the dose received by voxel i in the OAR at iteration k and of the desired dose for the same voxel,
wi=weights that determine how closely the minimization of "PHgr"( ) has to be carried out for voxel i in the PTV, and
xcex2i=weights that determine how closely the minimization of "THgr"( ) has to be carried out for voxel i in the OAR.
The dose received by a voxel i at iteration k is defined by       h    i          (      k      )        =            ∑      j        ⁢                  F        ij            ⁢              a        j                  (          k          )                    
wherein:
aj(k)=fluence for the jth pencil beam at iteration k and
Fij=element of the dose matrix F, giving the dose received in voxel i per unit fluence of pencil beam j.
An iteration consists of an attempt to minimize the cost function by increasing or decreasing the value of one randomly selected pencil beam fluence by a small amount.
Since the therapist only prescribes a range of maximum doses to be received by specific fractions of the OAR volume but cannot define which specific voxels are to receive a specific maximum dose, the desired doses si are calculated by a new procedure which takes into consideration the doses hj(k), for ixcex5S, after an optimization of the PTV only and the therapist""s specification.
The cost function has a second condition for inclusion of a particular OAR voxel i into the second summation. That condition is (hi(k)xe2x88x92si) greater than 0, which means that if a voxel receives a higher dose than desired in iteration k, it will be included in the summation. If a voxel already receives a lower dose than desired, that voxel does not enter into the computation of the cost function, or its derivatives, during iteration k. The number of terms in the cost function, then, changes dynamically and results in the xe2x80x9cadaptivexe2x80x9d nature of the new algorithm. The use of dynamic changes in the cost function of Simulated Annealing is another novelty of the ASA method.
An alternative method and apparatus for solving the numerical optimization problem comprises a computer running a form of the ASA iterative algorithm in which the analytic functions "PHgr"( ) and "PHgr"( ) are the quadratic forms
"PHgr"(hi(k),di)=(hi(k)xe2x88x92di)2
and
"THgr"(hi(k),si)=(hi(k)xe2x88x92di)2.
Accordingly, several objects and advantages of the invention are:
A primary object of the invention is to provide optimization of dose delivery in FMRT/FMRS procedures.
Another object of the invention is to provide radiation doses in the OAR that match or are less than the range of radiation doses specified by a therapist for the different voxels of the OAR.
Another object of the invention is to match radiation doses in the PTV to the desired dose in that volume as well as physical laws allow.
Another object of the invention is to provide optimization that allows the therapist to explore the space of available solutions in a clinically useful time.
An advantage of this invention is that it can provide optimizations for FMRT/S without over-constraining the PTV specifications.
Another advantage of this invention is that it can use analytic cost functions and respond to a detailed specification of which fractions of the OAR volume should receive no more than a certain dose range.
Yet another advantage of this invention is that optimizations can be obtained in relatively short computation times.
Other objects, advantages and novel features will be set forth in the following portions of the specification, wherein the detailed description is for the purpose of fully disclosing the preferred embodiment of the invention, without placing limitations thereon.